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ABSTRACT 


Waveform  modeling  of  radial  component  receiver  functions 
from  ANMO  (Albuquerque,  New  Mexico  Observatory)  for  three  source 
back  azimuths  (northwest,  southeast,  and  southwest)  has  been 
performed.  The  receiver  functions  were  derived  through  source 
equalization  deconvolution  of  merged  long  period  and  short  period 
digital  three  component  seismograms.  Derived  S-wave  velocity  models 
reflect  dominantly  intermediate  composition  granitic  rock  in  the 
upper  crust  (above  15  km  depth.  Vs  3.5  km/sec)  and  middle  crust 
(15-25  km  depth.  Vs  3. 5-3.7  km/sec).  Lower  crustal  shear 
velocities  of  approximately  3.75-3.85  km/sec  may  be  representative 
of  intermediate-to-maf ic  granulite  facies,  possibly  together  with 
previously  underplated  mafic  material  and  other  precursor  crustal 
rocks.  Shear  wave  attenuation  between  about  30-34  km  may  indicate  a 
lower  crustal  partial  melt  zone.  A  3-to-6  km  thick  interval  is 
interpreted  as  a  partial  melt  zone  in  the  upper  mantle  leading  into 
lessdepleted  spinel  peridotite  (Vs  =  4.25-4.35  km/sec)  near  37  km. 
Inversion  of  EPT-ALQ  interstation  dispersion  data  for  average  S-wave 
velocity  structure  produces  a  satisfactory  velocity  tie  to  the 
middle  and  lower  crust  portions  of  the  southwest  back  azimuth  model. 
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INTRODUCTION 


Receiver  function  analysis  has  become  a  widely  accepted 
method  of  modeling  shear  wave  velocity  structure  beneath 
isolated  digital  three-component  seismic  stations.  With 
broadband  data,  this  approach  offers  resolution  of  l-to-3  km 
thick  subsurface  layers  to  depths  near  60  km.  Furthermore, 
analyzing  off  azimuth  teleseismic  arrivals  recorded  at  a 
single  station  allows  for  modeling  of  structure  within  a 
maximum  radius  of  about  35  km  at  these  upper  mantle  depths. 
Therefore,  determination  of  detailed  crustal  and  upper  mantle 
structure  is  possible  through  receiver  function  studies. 

Previous  seismic  projects  conducted  in  the  Rio  Grande 
rift  have  provided  information  on  local  crustal  thickness, 
average  crustal  and  upper  mantle  shear  structure,  and  local 
P-wave  velocity  structure  within  the  southern  and  central 
rift  province.  Therefore,  in  the  first  part  of  this  thesis  a 
receiver  function  study  is  presented  in  which  shear  velocity 
structure  beneath  the  Albuquerque,  New  Mexico  seismic 
observatory  (ANMO)  is  modeled.  ANMO  is  located  on  the 
eastern  margin  of  the  Albuquerque  basin  in  the  Rio  Grande 
rift  (Figure  1).  The  objectives  of  this  study  are  to  (1) 
determine  whether  or  not  local  crustal  anomalies  exist,  (2) 
constrain  the  shear  structure  of  both  the  crust-mantle 
interface  and  the  upper  mantle,  and  (3)  propose  probable 
lower  crustal  and  upper  mantle  compositions  supported  by  the 
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Figure  1 . 


Rio  Grande  rift  map  depicting  the  Albuquerque 
basin  and  station  ANMO  (modified  from  Kelley, 
1979  )  . 
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geophysical  information. 

After  selecting  three-component  seismograms  of 
significant  teleseismic  events  (Mb  >=5.9)  digitally  recorded 
at  ANMO  during  1983  and  1984,  these  long  and  short  period 
data  are  merged  using  the  technique  of  Harvey  and  Choy 
(1981).  By  combining  the  long  and  short  period  spectra, 
records  of  broadband  displacement  (0.01  to  above  1.0  Hz 
bandwidth)  are  attained,  thus  improving  the  sensitivity  to  P- 
to-S  converted  phases  which  occur  within  the  crust  and  upper 
mantle.  Receiver  functions  are  then  derived  from  these 
records  using  the  processing  methods  of  Langston  (1979)  and 
Owens  (1984).  The  receiver  functions  are  then  modeled 
through  both  forward  raytracing  and  inversion  techniques  to 
attain  estimates  of  shear  velocity  structure  along 
propagation  paths. 

"^he  second  phase  of  this  thesis  is  a  reinvestigation  of 
the  Rayleigh  wave  dispersion  study  originally  conducted  by 
Sinno  and  Keller  (1986).  In  this  previous  study,  an  average 
S-wave  velocity  model  for  the  Rio  Grande  rift  axis  was 
derived  based  on  Rayleigh  wave  dispersion  along  the 
propagation  path  between  the  EPT  (El  Paso,  Texas)  and  the  ALQ 
(Albuquerque,  New  Mexico)  seismic  stations.  The  reanalysis 
of  data  used  in  the  original  study  involves  implementing 
phase  matched  filtering  to  improve  estimates  of  group  and 
phase  velocity,  and  time  variable  filtering  to  smooth  higher 
mode  and  multipath  interference  for  interstation  Green's 
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functions  determined  from  multi-event  data.  Simultaneous 
fits  of  dispersion  curves  derived  from  filtered  Green's 
functions  are  inverted  for  S-wave  velocity  structure  using 
the  inversion  scheme  of  Russell  (1980).  Both  biased  and 
unbiased  starting  velocity  models  for  the  inversion  are  used, 
where  the  biased  model  is  based  on  the  teleseismic  waveform 
modeling  results. 

This  study  presentation  is  organized  into  4  chapters. 
Chapter  1  covers  the  geologic  setting  of  the  Rio  Grande  rift, 
and  discusses  the  local  geology  for  the  Albuquerque  basin 
area.  Chapter  2  provides  a  summary  of  relevant  past 
geophysical  studies  in  the  Rio  Grande  rift.  A  discussion  of 
the  theory  and  methods  of  both  the  receiver  function  and 
surface  wave  studies  is  made  in  Chapter  3.  Lastly,  Chapter  4 
presents  the  processing,  modeling,  and  interpretation  results 
for  both  studies. 

A  summary  of  all  computer  code  used  in  this  project  is 
given  in  the  Appendix.  This  summary  includes  the  code  name, 
the  name  of  the  author,  modifications  made,  and  a  brief 
description  for  each  program. 


CHAPTER  1;  GEOLOGIC  SETTING 


RIO  GRANDE  RIFT 

The  Rio  Grande  rift  is  an  extensional  tectonic  feature 
that  extends  from  Leadville,  Colorado  to  Presidio,  Texas  and 
Chihuahua,  Mexico,  a  distance  of  more  than  1000  miles  (1600 
km)  (Figure  2)  (Baldridge  et  al . ,  1984).  Positioning  of  the 
rift  features  has  been  controlled  locally  by  underlying 
crustal  weaknesses  formed  during  the  many  orogenic  events 
that  occurred  from  the  late  Precambrian  to  the  early  Cenozoic 
(Chapin,  1979;  Ramberg  and  Smithson,  1975;  Chapin  and  Seager , 
1975).  Following  Laramide  tectonism,  the  compressional 
stress  field  ended  (Chapin,  1974)  and  subduction  of  the 
Farallon  Plate  beneath  the  North  American  Plate  resulted  in 
widespread  calc-alkalic  volcanism  (Cook  et  al .  ,  1979  ). 

Recent  studies  (Morgan  and  Golombek,  1984;  Seager  et 
al .  ,  1984  ;  Morgan  et  al .  1986)  demonstrate  two  phases  of 
extension  in  the  northern  and  southern  rift.  The  first  phase 
of  Cenozoic  extensional  strain  was  initiated  about  30  Ma . 
Crustal  sagging  along  the  developing  rift  formed  broad, 
shallow  basins  containing  mafic  flows  and  volcanic  ash  beds 
intercalated  with  alluvial  fill.  Complex  normal  faulting 
further  characterized  this  phase  (Chapin,  1979).  This 
transition  from  a  subduction  to  an  extensional  regime  is 
marked  by  a  change  from  calc-alkalic  to  basaltic  andesite 
magmatism.  This  phase  was  followed  by  a  magmatic  lull  which 
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Figure  2 


Tectonic  map  of  the  Rio  Grande  rift  (modified  from 
Baldridge  et  al . ,  1984). 
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lasted  approximately  between  26  Ma  to  13  Ma  in  the  southern 
rift  and  about  20  Ma  to  13  Ma  in  the  northern  rift  portion 
(Seager  and  Morgan,  1979). 

The  second  phase  of  extension  began  about  13  Ma .  At 
this  point  the  crust  became  critically  stretched,  allowing 
partial  melting  in  the  upper  mantle  bulge  beneath  the  rift. 
This  period  is  marked  by  emplacement  of  alkali  olivine 
basalts  (Seager  and  Morgan,  1979).  According  to  Chapin 
(1979),  mantle  upwelling  has  resulted  in  uplifts  of  about 
1100  m  since  the  middle  Miocene  in  the  Rocky  Mountains  and 
adjacent  areas.  Geophysical  and  regional  uplift  observations 
in  the  southern  Rio  Grande  rift  suggest  that  there  is  no 
direct  cause  and  effect  relationship  between  extension  and 
subsidence  as  would  be  expected  from  crustal  stretching 
models  with  a  constant  volume  crust,  so  that  underplating  of 
material  at  the  base  of  the  crust  is  required  to  explain  the 
regional  elevation  history  (Keller  et  al . ,  1989).  Bolson 
fill  thicknesses  of  1500  to  2500  m  are  common  in  basins  of 
the  southern  rift  (Seager  and  Morgan,  1979). 

Modern  activity  in  the  rift  is  evidenced  by  abundant 
fault  scarps  cutting  Pleistocene  deposits,  high  heat  flow, 
recent  elevation  changes,  recent  magma  bodies  and 
microseismicity  (Chapin,  1979). 

ALBUQUERQUE  BASIN  GEOLOGY 


Figure  1  depicts  the  location  of  the  Albuquerque  basin. 
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This  is  one  of  the  major  sedimentary  basins  along  the  length 
of  the  Rio  Grande  rift,  with  approximate  dimensions  of  100  mi 
(160  km)  north-south  and  40  mi  (64  km)  east-west.  The 
Precambrian  crystalline  basement  rocks  of  the  Sandia-Manzano- 
Los  Pinos  uplifts  bound  the  basin  on  the  eastern  margin.  The 
Albuquerque  basin  consists  of  two  sub-basins.  These  are  the 
North  and  South  Graben  Blocks  shown  in  Figure  3.  These 
basins  reach  maximum  depths  of  about  24,000  ft  (7317  m)  in 
the  North  Graben  and  23,000  ft  (7012  m)  in  the  South  Graben. 
This  basin  fill  is  Tertiary  continental  clastic  sedimentary- 
rocks,  overlying  up  to  11,000  ft  (3354  m)  of  Paleozoic  and 
Mesozoic  pre-rift  sedimentary  basement.  The  Precambrian 
crystalline  basement  consists  of  1300  Ma .  granitic  rocks  and 
1700  Ma .  metamorphic  rocks  (Condie,  1976;  Kelley,  1977; 

White,  1979)  (Russell  and  Snelson,  1990). 

CENTRAL  RIFT  MAGMATISM 

Figure  4  shows  the  main  volcanic  fields  of  the  central 
Rio  Grande  rift.  Lavas  from  the  northern  Albuquerque-Belen 
basin  consist  of  lower-alkali  olivine  tholeiite,  similar  to 
basalt  from  the  northern  rift.  In  the  central  Albuquerque- 
Belen  basin,  lavas  are  olivine  tholeiite,  basaltic  andesite, 
and  alkali  olivine  basalt  of  relatively  restricted 
compositional  range.  Experimental  data  show  that  the  olivine 
tholeiite  basalt  was  derived  from  partial  melting  (15-20%)  in 
the  upper  mantle,  and  that  the  alkali  olivine  basalt 


Figure  3.  E-W  geologic  cross  sections  of  l.he  Albuquerque 
basin  (from  Russell  and  Snelson,  1990). 
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originated  from  a  lesser  degree  of  partial  melting  at  50-70 
km  depth  (Baldridge,  1979  ).  Perry  et  al .  (1988)  suggest  that 
shifts  in  composition  from  alkali  basalts  to  tholeiites  was 
due  to  a  significant  shallowing  in  the  zone  of  magma  genesis 
following  a  thermal  conversion  of  lithospheric  mantle  to 
asthenosphere  between  8  and  4  Ma. 


CHAPTER  2:  PREVIOUS  GEOPHYSICS 


Gravity,  seismic  refraction,  seismic  reflection,  surface 
wave  dispersion,  teleseismic  travel-time  delay,  and  heat  flow 
studies  have  been  conducted  within  the  Rio  Grande  rift.  Such 
measurements  have  evidenced  crustal  thinning,  anomalously  low 
upper  mantle  density,  thick  accumulations  of  sediment  fill  in 
rift  basins,  and  seismic  anomalies  indicative  of  intracrustal 
magma  bodies. 

Rio  Grande  rift  gravity  modelling  efforts  have  revealed 
three  main  components  in  the  gravity  field  (Callender,  et  al. 
1989).  A  long  wavelength  gravity  low  is  interpreted  to  be 
due  to  lithospheric  thinning  (e.g.,  Bridwell,  1976;  Ander, 
1981).  An  axial  gravity  high  component  is  demonstrated  to  be 
due  to  crustal  thinning  (Daggett  et  al .  1986,  Decker  and 
Smithson,  1975;  Cordell,  1976,  1978;  Ramberg  et  al . ,  1978). 
This  anomaly  narrows  from  south  to  north,  indicating  rift 
extension  also  decreases  from  south  to  north  (Cordell,  1982). 
Lastly,  a  series  of  short  wavelength  gravity  lows  are  due  to 
thick  sediment  accumulations  in  rift  basins  (Callender,  et 
al.  198"). 

Seismic  refraction  studies  within  the  rift  (Table  1) 
include  the  works  of  Olsen  et  al .  (1979)  and  Sinno  et  al . 

(1986)  (published  results  of  Sinno,  1984).  The  1979  study 
involved  an  axial  profile  in  the  central  to  northern  rift 
segment.  The  1986  study  was  conducted  in  the  southern  Rio 
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Grande  rift  (Figure  5).  Olsen  eL  al  .  (  1979)  determined  an 

upper  crustal  P-wave  velocity  of  6.0  k.m/sec,  6.4  km/sec  for 
the  lower  crust,  and  an  upper  mantle  velocity  (Pn)  of  7.6 
km/sec,  with  a  crustal  thickness  of  34  km.  P-wave  velocities 
of  5.92  km/sec  for  the  upper  crust,  6.06  for  the  mid-crust, 
6.55  km/sec  for  the  lower  crust,  and  7.7  km/sec  (Pn)  for  the 
upper  mantle  were  published  by  Sinno  et  al .  (  1986  ).  Crustal 

thickness  was  determined  to  be  28  km  ne^r  the  southern  New 
Mexico  border  and  32  km  at  the  extreme  northern  end  of  the 
profile  (approximately  33.5  degrees  N  latitude).  This 
thinned  crust  distinguishes  the  rift  from  adjacent  the 
Colorado  Plateau,  Rocky  Mountains,  and  Great  Plains 
provinces.  Furthermore,  the  lower  observed  Pn  velocity 
distinguishes  the  rift  from  the  adjacent  Basin  and  Range 
province . 

Brown  et  al .  (1979,1980)  reported  results  of  seismic 

reflection  surveys  conducted  within  the  rift  by  the 
Consortium  for  Continental  Reflection  Profiling  (COCORP) 
(Figure  6).  A  strong  midcrustal  reflector  correlates  with 
the  tabular  magma  body  beneath  the  Socorro,  New  Mexico  region 
previously  detected  through  local  microearthquake  studies  by 
Sanford  et  al.  (1973  )  and  Brocher  (1981b)  (de  Voogd  et  al  .  , 
1986).  A  discontinuous  reflection  at  11-12  sec  agrees  with 
the  placement  of  the  M-discontinuity  derived  from  seismic 
refraction  by  Olsen  et  al .  (1979)  (Callendar  et  al.,  1989). 

The  crust-mantle  transition  zone  (between  33  to  38  km  depth, 
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Map  showing  previous  surface  wave  and  refraction 
work:  EPT-ALQ  surface  wave  path  (Sinno  and  Keller 
1986  study)  and  refraction  lines  (  a  =  Sinno  and 
Keller,  1986  axial  line;  b  =  Olsen  et  al . ,  1979 
axial  line)  (From  Sinno  and  Keller,  1986). 


Figure  5. 
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according  to  Toppozada  and  Sanford,  1976,  and  Olsen  et  al., 
1979)  was  revealed  to  be  a  complex  zone  varying  both 
laterally  and  vertically  on  a  scale  of  a  few  kilometers,  as 
evidenced  by  discontinuous  reflectors,  some  of  which  have  a 
layered  character  (Brown  et  al . ,  1979). 

Average  S-wave  velocity  models  for  propagation  paths 
along  the  rift  axis  and  across  the  cen'  ul  rift  have  been 
derived  from  Rayleigh  wave  studies.  These  studies  include 
Sinno  and  Keller  (1986),  and  Schlue  et  al .  (1986). 

Sinno  and  Keller  (1986),  using  seismic  records  from  5 
earthquakes  located  in  the  eastern  Pacific  between  1977  and 
1982,  analyzed  Rayleigh  wave  dispersion  between  the  EPT  (El 
Paso,  Texas)  and  ALQ  seismic  stations  (Figure  5).  Group 
dispersion  curves  were  determined  using  the  moving  window 
technique  (Landisman  et  al . ,  1969)  and  phase  dispersion 
curves  were  determined  using  the  cross-multiplication 
technique  (Bloch  and  Hales,  1968).  Group  and  phase 
velocities  were  inverted  simultaneously  using  a  first  order 
difference  inversion  scheme  (Russell,  1980;  Braile  and 
Keller,  1975)  to  determine  an  earth  model  which  fit  both  data 
sets.  Next,  P-wave  velocities  derived  from  seismic 
refraction  (Sinno  et  al.,  1986)  were  used  to  constrain  a 
surface  wave  velocity  inversion  using  the  spectral  whitening 
technique.  These  results  along  with  those  from  the 
refraction  stuf ^es  appear  in  Table  1.  Here,  an  average 
crustal  thickness  of  32  km  along  the  rift  axis  with  an 
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average  upper  mantle  S-wave  velocity  of  4.25  km/s  is 
indicated.  Poisson  ratios  calculated  based  on  these  velocity 
data  indicate  that  lithospheric  rigidity  in  the  southern  rift 
decreases  with  depth,  indicating  higher  than  normal  upper 
mantle  temperatures  (Sinno  and  Keller,  1986). 

Schlue  et  al .  (1986)  determined  an  upper  crustal 
S-wave  velocity  model  for  the  Albuquerque-Belen  Basin  for  the 
path  between  stations  RKO  and  TMG  shown  in  Figure  7.  Three 
Nevada  Test  Site  (NTS)  events  recorded  between  1983  and  1984 
were  used.  Fundamental  mode  Rayleigh  wave  phase  velocities 
were  obtained  and  inverted  to  yield  an  S-wave  velocity  model 
for  the  upper  20  km  of  the  basin.  This  model  contained  two 
low  velocity  zones  near  6  km  depth  and  18  km  depth.  The 
velocity  reversal  near  18  km  depth  was  considered  due  to  the 
thermal  effects  of  the  midcrustal  magma  body  previously 
identified  through  earthquake  studies  by  Sanford  et  al . 
(1973).  The  low-velocity  zone  detected  at  near  6  km  depth, 
based  on  high  heat  flow  values  of  110  mW/m  along  the  rift 
axis,  was  believed  to  be  possibly  due  to  the  effects  of  an 
upper  crustal  magma  body. 

Davis  et  al.  (1984)  analyzed  P-wave  travel-time  delay 
data  from  40  teleseismic  events  recorded  at  20  seismic 
stations  between  Moab,  Utah  and  Odessa,  Texas  between 
December  1982  and  January,  1983.  The  recording  station  array 
crossed  the  Rio  Grande  rift  axis  (south  of  Santal  Fe,  New 
Mexico)  at  a  45  degree  angle  along  an  azimuth  where 
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Figure  7.  Map  showing  RKO-TMG  surfcce  wave  path  (from 
Schlue  et  al . ,  1986). 
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teleseismic  arrivals  frequently  occur  (Figure  8).  The  study 
found  maximum  P-wave  delays  of  1.5  sec  for  the  central  rift. 
These  delays  were  interpreted  to  be  caused  by  an  upper  mantle 
low-velocity  zone  within  the  rift  boundary  in  the  depth  range 
70-100  km,  with  a  velocity  contrast  of  -8%  (Davis  et  al . , 

1984  )  . 

Estimates  of  heat  flow  from  a  deep  petroleum  test  in  the 
Albuquerque-Belen  basin  have  been  made  by  Reiter  et  al . 
(1986).  A  regional  heat  flow  of  —  77  mW/m^along  with  a  NE-SW 
trend  of  higher  heat  flow  (92-110  mW/m^ )  across  the  region 
were  determined.  The  background  heat  flow  is  consistent  with 
simple  models  of  crustal  thinning.  Higher  heat  flow 
components  are  believed  due  mainly  to  hydrothermal  fluids 
(largely  groundwaters  but  possibly  magma  in  some  areas) 
transporting  heat  upward  along  crustal  fracture  zones  (Reiter 
et  al . ,  1986).  Heat  flow  for  the  entire  rift  region  has  been 
mapped  by  Seager  and  Morgan  (1979)  (Figure  9). 
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Figure  8.  Portable  seismic  recording  station  locations  for 
teleseismic  P-delay  study  (from  Davis  et  al,, 
1984  )  . 
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Heat  flow  map  of  the  Rio  Grande  rift  region  (from 
Seager  and  Morgan,  1979,  and  Swanberg,  1979). 


Figure  9. 


CHAPTER  3:  THEORY  AND  METHODS 


This  thesis  involves  two  seismic  investigations  in  the 
Rio  Grande  rift  province:  a  receiver  function  study  using  the 
ANMO  seismic  station,  and  a  surface  wave  dispersion  study 
reanalysis  utilizing  data  from  the  EPT  and  ALQ  seismic 
observatories.  First,  the  methods  of  receiver  function 
analysis  as  applied  here  will  be  discussed.  Then  a 
discussion  of  the  concepts  behind  the  surface  wave  processing 
and  modeling  will  be  presented. 


RECEIVER  FUNCTION  STUDIES 
The  analysis  of  converted  phases  in  teleseismic 
waveforms  for  determination  of  local  earth  structure  beneath 
an  isolated  receiver  has  been  well  documented  {e.g.  Bath  and 
Stefansson,  1966;  Burdick  and  Langston,  1977;  Owens,  1984, 
1988).  Early  studies  attempted  to  resolve  earth  structure 
beneath  long-period  WWSSN  (Worldwide  Standardized  Seismograph 
Network)  stations.  The  limited  resolution  of  long-period 
data  allowed  only  for  the  determination  of  simple  crustal 
models.  Ov;ons  (1984  )  first  demonstrated  the  effectiveness  of 
using  broadband  digitally  recorded  data  to  resolve  detailed 
crustal  and  upper  mantle  structure. 

This  present  study  is  a  departure  from  prev'ious 
experiments  in  that  ANMO  lacks  broadband  recording 
instruments.  This  requires  beginning  with  the  digitally 
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recorded  three-component  long  and  short  period  seismograms 
from  ANMO,  and  merging  each  corresponding  long  and  short 
period  component  together  into  a  composite  record  with 
broadband  frequency  content.  Then,  proceeding  as  usual,  the 
merged  N-S  and  E-W  component  seismograms  are  rotated  with 
respect  to  their  source  back  azimuth,  and  source  equalization 
deconvolution  (Langston,  1979)  is  applied  to  isolate  receiver 
functions  from  the  records  of  total  horizontal  (radial  and 
tangential)  displacement.  These  receiver  functions  are  then 
gathered  into  stacking  suites  based  on  common  source  distance 
and  back  azimuth.  Each  suite  is  then  stacked  to  yield  a  mean 
receiver  function  for  the  particular  back  azimuth.  Synthetic 
waveform  matching  of  the  mean  radial  component  receiver 
functions  is  then  performed  in  both  forward  and  inverse 
modelling  schemes.  Langston’s  3-D  raytracing  method  (1979) 
is  implemented  in  both  approaches,  where  structure  beneath 
the  receiver  is  parameterized  with  a  stack  of  planar, 
homogenous  layers  of  finite  thickness  and  arbitrary  strike 
and  dip  (horizontal  layers  are  used  for  this  study).  The 
objective  is  to  reach  an  estimate  of  the  shear  velocity 
distribution  beneath  the  seismic  station  by  attaining  a  best- 
fit  waveform  match  between  the  observed  and  synthetic 
waveforms,  the  latter  being  the  result  of  raytracing  through 
the  parameterized  earth  model.  The  final  earth  model  yields 
the  best  waveform  match,  and  is  a  modified  version  of  a 
starting  model  into  which  all  available  a  priori  information 
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has  been  incorporated. 

DATA  SELECTION 

Following  the  criteria  outlined  by  Owens  (1984),  the 
following  requirements  are  met  by  the  chosen  data  set;  (1) 
Only  teleseismic  events  (distance  from  the  source  epicenter 
to  the  receiver  is  between  35-95  degrees)  are  used.  This 
guarantees  that  the  arriving  energy  will  impinge  steeply 
beneath  the  receiver  structure,  thereby  generating  an 
impulsive  vertical  response  within  the  structure. 

(2)  To  insure  significant  energy  reaches  the  receiver  for 
effective  modeling,  only  events  with  Mb  >=  5.9  are  used  for 
sources . 

PRELIMINARY  DATA  PROCESSING 

Prior  to  receiver  function  calculation,  initial 
processing  steps  must  be  performed  on  the  raw  digital  data. 
For  the  case  of  digitally  recorded  three-component  WWSSN  long 
and  short  period  data,  the  following  applies:  (1)  Using  the 
technique  of  Harvey  and  Choy  (1981),  the  long  and  short 
period  data  are  merged  together.  As  a  result,  the  merged 
records  have  a  bandwidth  from  0,01  Hz  to  above  1.0  Hz.  This 
significantly  yields  the  same  sensitivity  to  locally 
converted  P-to-S  phases  as  broadband  recording.  (2)  A 
rotation  is  applied  to  the  merged  horizontal  component  (N-S 
and  E-W)  records  in  order  to  determine  true  radial  and 
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tangential  receiver  displacement  with  respect  to  the 
particular  teleseismic  source  back  azimuth. 

Merging  the  digital  data  involves  first  aligning  the 
long  and  short  period  records  in  time.  Originally,  the  long 
period  seismogram  is  digitally  sampled  at  1  Hz,  and  the  short 
period  seismogram  at  20  Hz.  Before  further  processing,  the 
long  period  record  is  resampled  at  20  Hz.  Given  the 
difference  in  the  long  and  short  period  seismometer  responses 
(Figure  10),  these  effects  must  be  removed  prior  to  merging. 
Therefore,  the  long  and  short  period  seismograms  are 
separately  Fourier  transformed,  and  the  instrument  frequency 
responses  ar-  ■  ,-moved  in  a  divisional  deconvolution.  The 
effects  c*  _apid  instrument  gain  are  removed  before  merging 
by  an;;  lying  a  low-pass  filter  to  the  long  period  instrument 
cCi.rected  spectra  and  a  high-pass  filter  to  the  corrected 
short  period  spectra  (Figure  11).  An  inverse  Fourier 
transform  returns  the  long  and  short  period  time  series,  free 
of  instrument  effects.  These  are  then  summed  to  form  the 
resultant  broadband  displacement  seismograms  (Harvey  and 
Choy ,  1981 ) . 

Figure  12  shows  the  original  N-S  long  period  and  short 
period  records  of  a  particular  Honshu  teleseism  together  with 
the  corresponding  merged  broadband  record  used  in  this  study. 
The  inclusion  of  the  long  period  response  makes  possible  the 
extraction  of  velocity  information  at  lower  crustal  to 
subMoho  depths. 


AMPLITUDE  o  AMPLITUDE 


Long  period  SRO 


Figure  10.  SLaiion  ANMO  ( a )  short,  pe  riod,  and  (b)  long 

period  SRO  instrument,  responses  (from  Harvey  and 
Choy,  1981'!. 


Figure  11.  Illustration  of  frequency  cutoffs  for  low-pass 
(long  period)  and  high-pass  (short  period) 
filters  which  are  applied  to  ANMO  digital  data 
prior  to  merging  (from  Harvey  and  Choy,  1981). 


xample  of  data  merging.  LP  (long  period)  and  SP 
short  period)  displacement  records  are  merged  to 
reduce  the  BB  (broadband)  seismogram. 


30 


At  this  point  in  the  analysis,  the  data  set  consists  of 
three-component  (N-S,  E-W,  and  Z)  broadband  digital  seismic 
records.  The  next  step  involves  rotating  the  horizontal 
component  seismograms  with  respect  to  the  particular  source 
back  azimuth.  The  rotated  output  consists  of  true  radial  and 
tangential  displacement  seismograms  for  a  propagation 
direction.  Once  rotated,  source  equalization  can  be 
performed  on  the  three-component  data  to  determine  the 
receiver  functions. 

RECEIVER  FUNCTION  CALCULATION 

Following  the  approach  used  by  Langston  (1979),  the 
total  seismograms,  D(t),  recorded  at  the  three-component 
digital  station,  can  be  expressed  in  the  time  domain  as 
convolutions  of  the  impulse  response  of  the  recording 
instrument,  I(t^  the  source  function,  S(t),  and  the  response 


of  the  local  earth  structure. 

E  ( t )  : 

Dv{t)  =  I(t) 

*  S(t) 

*  Ev ( t ) , 

Dr(t)  =  Kt) 

*  S(t) 

*  Er(t) , 

Dt(t)  =  I(t) 

*  S(t) 

*  Et(t) , 

where  v,  r,  and  t  subscripts  denote  vertical,  radial,  and 
tangential  components  of  motion,  respectively.  Determination 
of  Lhe  receiver  function  requires  the  removal  of  the  unwanted 
source  and  instrument  effects  to  reveal  the  response  of 
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structure  beneath  the  receiver, 
responses  are  removed  during  the 
the  broadband  total  displacement 
as  : 


Given  that  instrument 
data  merging  process,  then 
seismograms  are  expressed 


Dv( t )  =  S(t)  *  Ev( t) , 
Dr(t)  =  S(t)  *  Er(t) , 
Dt{t)  =  S(t)  *  Et(t) , 


with  the  source  function  and  earth  response  denoted  as  usual . 
Further  development  here  will  also  exclude  the  instrument 
response  term. 

Assuming  that  the  vertical  component  earth  response  for 
teleseismic  events  will  be  dominated  by  mainly  a  large  direct 
P  arrival  followed  only  by  minor  crustal  reverberations  and 
phase  conversions  (Langston,  1979),  the  earth  response  for 
the  vertical  component  of  motion  can  be  approximated  by  the 
Dirac  delta  function: 


Ev  { t )  ( t )  . 

This  implies  that  the  total  vertical  component  seismogram  can 
be  represented  as 


Dv { t )  ^  S ( t )  , 
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meaning  that  the  source  effects  which  should  be  removed  are 
contained  here.  This  allows  for  the  source  equalization 
deconvolution  scheme  developed  by  Langston  (1979),  in  which  a 
frequency  domain  deconvolution  is  applied  to  reveal  the 
radial  and  tangential  component  receiver  functions: 

Er(w)  =  Dr(w)  /  (S(w))  c:  Dr(w)  /  Dv(w), 

Et(w)  =  Dt(w)  /  {S(w))  —  Dt(w)  /  Dv{w). 

While  stable  for  noise-free  data,  this  divisional 
deconvolution  becomes  unstable  when  applied  to  band-limited 
seismograms  containing  noise  (Owens,  1984).  Therefore,  a 
minimum  allowable  amplitude  level  (waterlevel)  is  applied  to 
fill  spectral  troughs.  Lastly,  smoothing  with  a  Gaussian 
function  is  performed.  The  final  frequency  domain  expression 
for  the  radial  receiver  function  is 

Er(w)  =  [ Dr ( w) Dv ( w) /phi ( w ) ]G ( w ) , 


where 


phi(w)  =  max { Dv ( w) Dv ( w) ,  c -max[ Dv(w) Dv(w)  ] } ; 
c  =  waterlevel  parameter; 

G(w)  =  exp(-w  /4a  )=  Gaussian  smoother. 


Inverse  transformation  of  the  radial  and  tangential  sp^-'Ctra 
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yield  the  time  domain  receiver  functions,  Er(t)  and  Et ( t ) . 

In  the  above  expression  for  G{w) ,  parameter  a  controls 
the  width  of  the  Gaussian  pulse,  in  turn  controlling  the 
degree  to  which  high  frequency  noise  is  reduced.  For  the 
frequency  content  of  broadband  receiver  functions,  a  =  5  is 
an  appropriate  selection  (Owens,  1984).  Including  a 
waterlevel  parameter  (c  >  0)  reduces  noise  in  the  receiver 
function  due  to  spectral  troughs  and  high  frequency  noise 
present  in  the  vertical  componei'.t,  but  a  reduction  in 
relative  arrival  time  resolution  can  be  expected  with 
increased  c  magnitude.  Furthermore,  increased  c  values  lead 
to  truncation  of  higher  frequencies  in  the  vertical  component 
amplitude  spectrum,  which  can  introduce  sidelobes  into  the 
receiver  functions.  An  optimum  tradeoff  between  these 
effects  must  be  made,  which  the  processor  must  determine  by 
trial  and  error.  Appropriate  values  of  c  fall  in  the  range 
between  0.0001  to  1.0.  The  deconvolution  should  be  performed 
over  this  range,  and  the  optimum  result  chosen  (Owens,  1984). 

STACKING  RECEIVER  FUNCTIONS 

Quality  receiver  functions  are  gathered  into  stac)<.ing 
suites  based  on  relative  origin  distance  and  back  azimuth  and 
stacked  to  increase  signal  to  noise  ratio,  and  to  provide  an 
average  estimate  of  crustal  and  upper  mantle  response  along  a 
particular  propagation  path.  The  range  in  distance  variation 
within  a  single  stacking  suite  should  be  less  than  15  degrees 
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for  events  with  distance  >=  70  degrees,  and  less  than  10 
degrees  for  closer  events.  Back  azimuthal  variation  should 
be  kept  to  less  than  20  degrees.  Determination  of  best 
estimate  receiver  functions  for  a  range  of  back  azimuths 
(controlled  by  the  natural  distribution  of  teleseismic 
sources,  mainly  NW,  SE,  and  SW  back  azimuths  for  this  study) 
allows  for  modeling  of  horizontal  variations  in  shear  wave 
velocity  structure  beneath  the  station  (Owens,  1984). 

MODELING  RECEIVER  FUNCTIONS 

The  3-D  raytracing  method  of  Langston  (1977)  requires 
the  assumption  that  earth  structure  beneath  the  three- 
component  receiver  can  be  modeled  as  a  stack  of  homogeneous, 
planar  layers  of  varying  thickness  and  dip  orientation 
overlying  a  half -space.  Furthermore,  each  layer  is 
parameterized  with  an  S-wave  velocity,  a  P-wave  velocity,  and 
a  density.  Since  teleseismic  sources  are  used,  plane  waves 
are  assumed. 

P-waves  arriving  from  teleseismic  sources  impinge  at 
steep  angles  of  incidence  beneath  the  layered  stack  (Figure 
13).  As  the  energy  propagates  upward  through  the  layers,  P- 
to-S  conversions  occur  due  to  velocity  contrasts  at  layer 
interfaces.  These  shear  waves  are  recorded  on  the  horizontal 
component  seismograms  within  50  seconds  following  the  arrival 
of  the  direct  P  phase  for  interfaces  up  to  80  km  depth 
(Owens,  1984).  Langston's  method  assumes  all  multiples  are 
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between  the  free  surface  and  layer  interfaces  within  the 
stack;  hence,  no  intralayer  multiples  can  be  included.  The 
notation  of  Bath  and  Stefansson  (1966)  is  used  to  name 
specific  phases  produced:  (1)  Incident  waves  at  the  bottom  of 
a  layered  stack  are  represented  by  capital  letters  (P  or  S). 
(2)  Upgoing  waves  in  the  layered  stack  are  represented  by 
lower  case  letters  (p  or  s).  (3)  Downgoing  waves  in  the 
layers  (reflections  from  the  free  surface  )  are  represented 
by  capital  letters. 

(4)  Reflected  waves  off  the  top  of  an  interface  are 
represented  by  subscripts  indicating  the  number  of  the 
reflecting  interface  or  an  m  for  the  Moho  (Owens,  1984). 
Extensive  forward  modeling  was  used  in  this  study,  along  with 
implementation  of  an  inverse  scheme  developed  by  Owens. 

In  the  time  domain  inversion,  the  starting  model  is 
assumed  to  be  close  to  the  true  model,  and  making  a  Taylor 
series  expansion  about  the  initial  model  yields  the 
linearized  equations 

N 

dRi  =  Ri(obs)  -  Ri(r)k)  =  ^ 

kViSp*' 

Here,  subscript  i  denotes  the  ith  sample  of  the  observed  and 
synthetic  radial  receiver  functions,  Ri(ob5)  and  Ri(pk), 
where  pk  is  the  kth  initial  model  parameter,  and  dRi  is  the 
residual  between  the  observed  and  synthetic  waveforms.  The 
residual  is  minimized  using  the  L2  norm,  and  the  pk  are 
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iteratively  corrected  by  dpk  until  a  satisfactory  waveform 
match  is  achieved. 

Synthetics  are  calculated  using  Langston's  raytracing 
method,  as  in  the  forward  problem,  and  only  S-wave  velocity 
is  inverted  for.  The  inverse  operator  is  determined  using 
singular  value  decomposition,  and  the  partial  derivatives  are 
determined  by  a  finite  difference  scheme  (Owens,  1984). 

It  should  be  noted  here  that  both  forward  and  inverse 
modeling  efforts  are  exploited  in  this  study.  Initially, 
extensive  forward  modeling  offers  insight  into  relating 
observed  phases  with  specific  structural  zones  at  depth. 
Basing  the  intial  models  mainly  on  previous  seismic 
refraction  and  surface  wave  dispersion  results,  forward 
modeling  to  refine  the  models,  and  inverting  the  receiver 
functions  with  smoothed  versions  of  these  well  constrained 
starting  models  (implying  that  extreme  bias  is  removed  from 
the  starting  models  before  inversion)  increases  confidence  in 
the  final  derived  estimates  of  S-wave  velocity  structure. 


SURFACE  WAVE  ANALYSIS 

Rayleigh  wave  dispersion  in  the  Rio  Grande  rift  has  been 
previously  analyzed  by  Sinno  and  Keller  (1986)  .  They 

employed  the  two-station  method,  using  events  recorded  at 
both  the  EPT  and  ALQ  WWSSN  seismic  stations,  determining 
interstation  group  and  phase  velocities  from  fundamental  mode 
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Rayleigh  wave  data.  The  group  and  phase  velocities  were 
simultaneously  inverted  for  gross  shear  velocity  structure 
along  the  interstation  path. 

In  this  thesis,  the  data  used  in  the  1986  original  study 
is  reprocessed  with  a  scheme  which  involves  filtering 
seismograms,  computing  the  interstation  Green's  function,  and 
determining  the  self-consistent  phase  and  group  dispersion 
with  standard  errors  (Dean  and  Keller,  to  be  published).  The 
first-order  difference  inversion  scheme  (Russell,  1980)  is 
used  to  derive  shear  velocity  models  for  the  interstation 
path . 

DATA  SELECTION 

The  two  station  method  requires  that  in  order  to 
determine  interstation  shear  structure  via  Rayleigh  wave 
dispersion  analysis,  the  earthquake  sources  must  occur  within 
a  few  degrees  of  the  great  circle  path  between  the  stations 
(Knopoff  et  al . ,  1967).  All  5  events  used  in  the  1986  study 
met  this  requirement.  For  the  present  reanalysis,  paper 
seismogram  records  from  both  the  EPT  and  ALQ  stations  were 
obtained  and  the  surface  wave  trains  were  redigitized.  This 
data  set  consisted  of  records  for  3  of  the  original  5  events 
used  . 

SURFACE  WAVE  DATA  PROCESSING 


The  processing  scheme  devised  by  Dr.  E.  A.  Dean  (see 
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Appendix)  Lo  derive  interstaLion  group  and  phase  velocities 
from  multi-event  data  sets  was  implemented.  As  outlined  by 
Dean  and  Keller  (to  be  published),  this  scheme  involves  the 
following  steps  for  a  single  event: 

(1)  Resampling  the  raw  digitized  surface  wave  signal, 
interpolating  at  1  sample/sec, 

(2)  Windowing  the  relevant  portion  of  the  seismograms  is 
performed  to  reduce  the  number  of  samples  processed.  By 
default,  this  is  the  portion  with  apparent  velocities  between 
1.8  )cm/sec  and  4.6  km/sec. 

(3)  Optional  application  of  a  zero-phase  quasi-pin)c  filter 
(QPF)  which  smooths  the  spectrum  towards  constant  power  per 
octave . 

(4)  Application  of  the  multiple  filter  technique  (MFT)  to 
determine  initial  group  velocity  estimates.  A  matrix  of 
power  vs.  period  and  group  velocity  (arrival  time)  is 
displayed  on  the  computer  screen,  and  the  user  may  edit  the 
computer  selected  ridge  of  the  contoured  matrix  by  choosing  a 
different  local  maximum  or  deleting  a  selection. 

(5)  A  least  squares  polynomial  in  log  T  smooths  the  group 
velocity  function  determined  by  the  MFT. 
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(6)  Phase  matching  iteration  (PMI)  refines  either  phase  or 
group  velocity.  The  phase  phi(f)  is  determined  by  numerical 
integration  of  the  group  arrival  time  t(f): 

pKi  (f)  =  j -t  (f 

o 

A  "pseudo-autocorrelation  function"  (PAF)  is  derived  by  cross 
correlating  a  unit-amplitude  signal  of  phase  phi(f)  with  the 
seismogram.  This  PAF  is  windowed  in  time  and  its  phase 
spectrum  phi ' ( f )  is  used  to  improve  the  group  arrival  times 
by  the  inverse  of  the  above  equation, 

t' (f )  =  dphi’ (f )/df . 

This  technique  is  iteratively  applied  until  the  PAF  is 
symmetric . 

(7)  A  time  variable  filter  (TVF)  based  on  the  refined  group 
velocity  is  applied.  This  step  attenuates  higher  modes  and 
multipath  arrivals. 

(8)  The  filtered  signals  x(t)  (near  station)  and  y(t)  (far 
station)  are  then  used  to  determine  an  interstation  Green's 
function.  The  frequency  domain  convolutional  representation 
of  the  far  station  signal  can  be  expressed  as 
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Y(f }  =  G(f )  .  X(f )  , 

where  G(f)  is  the  spectrum  of  the  interstation  Green's 
function.  G(f)  can  be  recovered  by  the  spectral  ratio  method 
as  follows: 


G(f )  =  C(f )  /  A(f ) , 

where  C(f)  is  the  transformed  cross  correlation  of  x(t)  and 
y(t),  and  A(f)  is  the  transformed  autocorrelation  of  x(t). 
Spectral  holes  in  the  near  station  signal's  spectrum  can 
cause  instability  in  the  deconvolution.  The  deconvolution 
performed  as 


<G'  (f  )>  =  <C(f  )/[  lx(f  )  1  I  Y(f  )  1  ]^> 

decreases  this  effect,  where  <  >  on  the  right  side  of  the 
equation  means  windowing  the  PAF  of  g'(t),  and  ]x(f)|  and 
lY(f)|  are  the  amplitudes  of  the  filtered  seismogram  spectra. 
Inverse  transformation  of  G'(f)  yields  the  interstation 
Green's  function,  g'{t),  which  is  then  filtered  as  outlined 
in  the  steps  1-7.  Further  analysis  to  determine  the 
interstation  group  and  phase  velocities  is  as  follows: 

(9)  MFT  is  applied  to  the  filtered  Green's  function  to 
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determine  the  initial  estimate  of  group  velocity  (u(T)). 

(10)  A  least  squares  polynomial  of  variable  order  smooths  the 
MFT  estimated  group  velocity.  The  F-test  (Dean,  1986)  is 
used  to  determine  the  highest  order  significant  fit. 

Standard  errors  (/\u)  for  the  smoothed  velocities  are 
determined  from  the  covariance  matrix  of  the  polynomial  fit 
to  the  group  velocity  data. 

(11)  PMI  is  applied  to  both  fast  (u(T)  +  /\u )  and  slow  (u(t) 
-/\u)  group  velocities  and  the  two  refined  phase  velocities 
are  meaned  to  obtain  "raw  phase  velocity". 

(12)  A  simultaneous  least-squares  fit  to  the  raw  group  and 
phase  velocities  is  applied  as  described  by  Dean  (1986).  A 
modified  version  of  the  dispersion  relation 

u  =  v/[ l+(T/v) (dv/dT)  ] 

is  used  which  handles  the  non-linear  least  squares  analysis 
and  log-T  dependence  in  this  approach.  Wor)cing  with  group 
and  phase  slownesses  (y(x)  =  1/u,  z  =  1/v)  and  a  scaled 
expression  for  period  (x  =  A  -  3(logT),  with  A  and  B  such 
that  X  =  -1  for  the  minimumi  period  of  analysis,  and  x  =  1  for 
the  m.aximum  period  of  analysis)  yields 
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y(x)  =  z(x)  -  B(dz/dx). 

Phase  slowness  is  approximated  by 

z  {  X  )  =  ^ci  'Pi  (  X  )  , 

where  Pi(x)  is  an  ith  degree  polynomial  of  x. 
is 


Group  slowness 


y(x)  =  ^ci.Qi(x), 

where  Qi(x)  =  Pi(x)  -  B ( dPi ( x ) /dx ) . 

The  expansion  coefficients  are  solved  by  the  least  squares 
method  using  singular  value  decomposition,  and  the  order  of 
fit  is  determined  by  the  F-test.  Errors  (/\u,  ^v )  for  the 
fitted  group  (u)  and  phase  (v)  velocities  are  determined  from 
the  covariance  matrix  for  the  singular  value  decomposition  of 
the  design  matrix  (Dean  and  Keller,  to  be  published). 

MODELING  GROUP  AND  PHASE  VELOCITY  DISPERSION 

The  first-order  difference  inverse  technique  (Russell, 
1980)  is  used  to  derive  interstation  shear  velocity  models 
from  the  phase  and  group  velocities.  This  method 
simultaneously  minimizes  the  standard  least  squares  error  and 
the  magnitude  of  the  difference  between  adjacent 


solution 
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elements  (Twomey,  1977).  Expressed  mathematically, 

N  ,  £:t  i 

(i-A)Xek  +  A^(Xk-i  -Xk) 

krl 

is  minimized,  where  ek  is  the  error  between  observed  and 
predicted  values,;  xk  is  the  kth  element  of  the  solution 
(model)  vector,  x;  A  is  an  arbitrary  scale  between  0  and  1;  N 
is  the  number  of  observations;  M  is  the  number  of  model 
parameters.  In  the  above  equation,  the  second  term  acts  to 
smooth  the  optimum  solution  (Sinno  and  Keller,  1986). 
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CHAPTER  4:  SEISMIC  STUDIES  AND  RESULTS 

ANMO  RECEIVER  FUNCTION  STUDY 

In  this  section,  both  the  processing  and  velocity 
modeling  results  for  the  receiver  function  analysis  are 
presented.  Important  features  in  both  waveforms  and  velocity 
models  will  be  noted.  Geologic  interpretations  will  be 
addressed  in  the  following  section. 

Table  2  lists  the  teleseismic  events  used  as  sources  in 
this  study.  These  were  all  recorded  at  ANMO  between  1983- 
1984.  Based  on  their  spatial  distribution,  local  structure 
is  modeled  along  three  constrained  back  azimuths:  Honshu, 
Chile,  and  Tonga  events  provided  control  along  the  northwest, 
southeast,  and  southwest  back  azimuths,  respectively. 

Only  the  radial  component  receiver  functions  are  modeled 
in  this  study.  Although  tangential  motion  can  provide 
information  on  lateral  structural  variation  (e.g.,  dipping 
interfaces),  the  merged  tangential  receiver  functions 
demonstrate  no  characteristic  polarity  trends  with  azimuthal 
variation.  Therefore,  horizontal  layer  parameterizations  are 
made,  and  only  observed  radial  displacement  is  modeled. 

Figure  14  shows  the  northwest  back  azimuth  stacking 
suite  receiver  functions.  The  waterlevel  (c)  values  used  in 
the  deconvolution  process  are  listed  for  each  trace.  The 
resultant  mean  stacked  receiver  function  ir  given  at  the 
bottom  of  Figure  14.  Figure  15  is  the  mean  trace  again,  with 
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TABLE  2:  ANMO 

Teleseismic  Events 

Date 

Distance 

Back  Azimuth 

Depth 

Mb 

<degrees> 

<degrees> 

<km> 

Northwest 

Back  Azimuth 

01/01/84 

89 . 6 

311.8 

368. 

6.5 

03/06/84 

91.2 

307.5 

457  . 

6.2 

09/18/84 

86 . 5 

309.5 

48. 

6.6 

05/26/83 

83.7 

315.7 

24. 

6.8 

06/09/83 

83.9 

315.5 

31 . 

6.3 

08/25/83 

92.9 

314.8 

126. 

6.1 

Southwest 

Back  Azimuth 

06/01/83 

82.7 

243.5 

180. 

6.2 

09/17/83 

84.6 

245.6 

33. 

6.1 

08/30/83 

80.7 

242.2 

39. 

6.0 

10/17/83 

84 . 5 

240.1 

30. 

6.0 

12/03/83 

80.6 

243.6 

82. 

6.1 

06/15/84 

82.1 

244.6 

247  . 

6.1 

Southeast 

Back  Azimuth 

10/04/83 

69.9 

146.0 

15. 

6.4 

10/09/83 

69.6 

145.8 

16. 

5.9 

10/20/84 

69.7 

141.9 

213. 

5.9 

Receiver  functions  comprising  northwest  back 
azimuth  data  set.  A  =  source  distance;  BAZ  = 
back  azimuth;  c  =  waterlevel  value  used;  stack  = 
mean  stacked  trace. 


Figure  14. 
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c  t>  :> 


Figure  15. 


(a)  Mean  stacked  receiver  function  for  northwest 
back  azimuth,  and  (b)  +/-  1  standard  deviation 
data  bounds . 
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its  data  bounds.  The  data  bounds  are  determined  from  the 
variance  in  the  seismograms  within  the  stacking  suite,  and 
are  plotted  as  plus-or-minus  1  standard  deviation  about  the 
mean.  This  provides  a  statistical  check  of  the  waveform 
matching  in  that  a  satisfactory  match  between  the  synthetic 
and  observed  radial  receiver  functions  is  achieved  if  the 
synthetic  fits  within  these  bounds.  Since  traces  are 
normalized  to  the  amplitude  of  the  direct  P  arrival  before 
stacking,  a  zero  variance  results  about  this  first  arrival  in 
the  stacked  receiver  function.  The  nonzero  variance  at  the  P 
arrival  in  the  northwest  mean  trace  is  due  to  the  fact  that 
more  relative  weight  was  applied  during  stacking  to  the 
cleaner  traces  within  this  stacking  suite.  In  other  words, 
once  normalized,  the  cleaner  traces  were  given  a  greater 
relative  scaling  in  amplitude  than  the  noiser  seismograms, 
and  the  statistics  calculated  during  stacking  reflect  this. 

In  retrospect,  for  statistical  purposes,  it  would  have  been 
more  appropriate  to  weight  these  additively  rather  than 
multiplicatively .  However,  it  will  be  demonstrated  that  the 
followed  procedure  did  not  distort  these  results. 

Figures  16  and  17  show  the  southeast  back  azimuth 
stacking  suites  and  the  mean  stacked  trace  with  data  bounds, 
respectively.  The  similar  data  for  the  southwest  back 
azimuth  are  given  in  Figures  18  and  19.  Phases  seen  between 
2-5  seconds  (following  the  direct  P  arrival  at  0  seconds)  in 
these  mean  receiver  functionc  ^'Figures  17  and  19)  exhibit  a 
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Figure  16.  Receiver  functions  comprising  southeast  back 

azimuth  data  set,  A  =  source  distance;  BAZ  = 
back  azimuth;  c  =  water  level  value  used;  stack  = 
mean  stacked  trace. 


Figure  17.  (a)  Mean  stacked  receiver  function  for  southeast 

back  azimuth,  and  (b)  +/-  1  standard  deviation 
data  bounds . 
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Figure  18.  Receiver  functions  comprising  southwest  back 

azimuth  data  set.  A  =  source  distance;  BAZ  = 
back  azimuth;  c  =  waterlevel  value  used;  stack  = 
mean  slacked  trace. 


Figure  19.  (a)  Mean  stacked  receiver  function  for  southwest 

back  azimuth,  and  (b)  +/-  1  standard  deviation 

data  bounds . 
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similar  character.  They  are  significantly  different  from  the 
northwest  back  azimuth  receiver  function  (Figure  15)  within 
the  same  time  window.  In  each  mean  trace,  converted  phases 
arriving  between  13-20  seconds  are  multiples  from  the  crust- 
mantle  transition  zone.  Phases  following  these  arrivals  are 
conversions  from  deeper  structure  within  the  upper  mantle. 

The  final  radial  waveform  fit  for  the  northwest  back 
azimuth  is  given  in  Figure  20.  All  phases  between  0-27 
seconds  are  satisfactorily  matched,  as  indicated  by  both 
Figure  20-a  and  Figure  20-b.  The  corresponding  derived  S- 
wave  velocity  model  is  shown  in  Figure  21  and  Table  3.  The 
upper  crustal  portion  of  the  derived  model  indicates  a 
positive  shear  velocity  gradient  down  to  15  km  depth.  A 
slight  positive  velocity  step  near  28  km  depth  tends  to 
produce  the  Ps  phase  arriving  at  3.2  seconds  in  this  receiver 
function,  which  is  not  seen  in  the  other  two  receiver 
functions.  This  might  suggest  a  more  defined  mid-to-lower 
crust  division  toward  the  northern  rift  portion,  but  more 
data  would  be  needed  to  verify  this.  Just  beneath  this 
feature  at  about  30  km,  a  2  to  4  km  thick  zone  appears  which 
shows  a  slight  flattening  in  S-wave  velocity.  Upper  mantle 
velocities  are  encountered  near  37  km  depth,  and,  given  the 
+/-  2  standard  errors  in  the  velocity  estimates,  the  model 
indicates  the  transition  from  lower  crust  to  upper  mantle 
occurs  over  an  interval  of  approximately  3  to  6  km.  A 
positive  shear  wave  velocity  gradient  of  about  .01  km/sec/km 
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Figure  20.  Final  waveform  fit  for  northwest  back  azimuth. 

(a)  Dotted  trace  =  mean  receiver  function,  solid 
trace  =  derived  synthetic.  (b)  Dotted  traces  = 
upper  and  lower  data  bounds,  solid  trace  = 
derived  synthetic. 
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Figure  21.  Final  S-wave  velocity  model  for  northwest  back 
azimuth.  Error  bars  represent  +/-  2  standard 
errors . 
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TABLE  3;  Final  S-wave  Velocity  Model —  ANMO  NW  Back  Azimuth 


Layer  # 

Thickness 

<km> 

NW  S-wave  Velocity 
<km/sec> 

1 

1,75 

2.25 

2 

1.60 

3.50 

3 

1.50 

3.42 

4 

2.10 

3.57 

5 

2.00 

3.58 

6 

2.25 

3.61 

7 

2.25 

3.67 

8 

2.25 

3.45 

9 

2.00 

3.46 

10 

2.75 

3.70 

11 

2.25 

3.59 

12 

2.90 

3.60 

13 

1.40 

3.56 

14 

2.00 

3.87 

15 

2.50 

3.85 

16 

2,50 

3.77 

17 

2.80 

3.77 

18 

3.00 

4.11 

19 

3.00 

4.40 

20 

2.00 

4.08 

21 

2.00 

4.32 

22 

2.00 

4.37 

23 

2.00 

4.36 

24 

2.50 

4.43 

25 

2.50 

4.27 

26 

2.50 

4.37 

27 

2.50 

4.42 

28 

2.50 

4.43 

29 

0.00 

4.44 
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beneath  the  crust  mantle  transition  continues  to  over  60  km 
depth  in  this  model. 

Figures  22,  23,  and  Table  4  give  the  final  waveform 
modeling  results  and  final  velocity  model,  respectively,  for 
the  southeast  back  azimuth.  This  model  has  a  smooth  positive 
velocity  gradient  to  a  depth  near  28  km.  As  in  the  northwest 
model ,  shear  velocities  seem  to  be  attenuated  close  to  30  km 
depth,  but  when  the  standard  errors  in  the  southeast  velocity 
model  at  this  depth  are  considered,  verification  of  this 
particular  feature  is  difficult.  A  gradational  Moho  shear 
boundary  near  36  km  depth  determined  in  the  inversion  yields 
a  significantly  good  waveform  fit  for  this  back  azimuth. 

Here  again,  shear  velocities  within  the  upper  mantle  increase 
with  depth. 

Lastly,  Figure  24  is  the  derived  synthetic  fit  for  the 
southwest  back  azimuth.  Figure  25  and  Table  5  give  the 
derived  S-wave  velocity  model.  Upper  crustal  velocities 
increase  to  a  depth  of  15  km.  Interpretation  of  the  slight 
velocity  reversal  between  15  and  20  km  may  not  be  justifiable 
due  to  the  overlap  in  error  bars.  Again,  flattening  to 
reversing  of  shear  velocity  is  noted  at  30  to  32  km  depth.  A 
transitional  shear  boundary  grading  int  upper  mantle  near  36 
km  depth  produces  a  statistically  acceptable  match  for  the 
crust  mantle  multiples.  As  in  the  other  models,  the  crust 
mantle  transition  occurs  over  an  interval  of  3  to  6  kmi . 
Although  this  model  shows  a  strong  positive  velocity  step 


Figure  22. 


Final  waveform  fit  for  southeast  back  azimuth. 

(a)  Dotted  trace  =  mean  receiver  function,  solid 
trace  =  derived  synthetic.  (b)  Dotted  traces  = 
upper  and  lower  data  bounds,  solid  trace  = 
derived  synthetic. 


DEIRTM  <  U.  m  > 


Figure  23  .  Final  S-wave  velocity  model  for  southeast  back 
azimuth.  Error  bars  represent  +/-  2  standard 


errors 
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TABLE  4:  Final  S-wave  Velocity  Model —  ANMO  SE  Back  Azimuth 


Layer  # 

Thickness 

<km> 

SE  S-wave  Velocity 
<km/sec> 

1 

0.80 

2.51 

2 

1.60 

2.61 

3 

1.50 

3.56 

4 

2.10 

3.56 

5 

2.00 

3.20 

6 

2.25 

3.58 

7 

2.25 

3.32 

8 

2.25 

3.55 

9 

2.00 

3.52 

10 

2.75 

3.59 

11 

2.25 

3.59 

12 

2.90 

3.67 

13 

1.40 

3.80 

14 

2.00 

3.66 

15 

2.50 

3.58 

16 

2.50 

3.57 

17 

2.80 

3.89 

18 

3.00 

4.29 

19 

3.00 

4.28 

20 

2.00 

4.27 

21 

2.00 

4.22 

22 

2.00 

4.36 

23 

2.00 

4.25 

24 

2.50 

4.36 

25 

2.50 

4.21 

26 

2.50 

4.49 

27 

2.50 

4  .  36 

28 

2.50 

4  .  54 

29 

0.00 

4 .44 

Figure  24.  Final  waveform  fit  for  southwest  back  azimuth. 

(a)  Dotted  trace  =  mean  receiver  function,  solid 
trace  =  derived  synthetic.  (b)  Dotted  traces  = 
upper  and  lower  data  bounds,  solid  trace  = 
derived  synthetic. 
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Figure  25 .  Final  S-wave  velocity  model  for  southwest  back 
azimuth.  Error  bars  represent  +/-  2  standard 


errors  . 
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TABLE  5;  Final  S-wave  Velocity  Model —  ANMO  SW  Back  Azimuth 


Layer  # 

Thickness 

<km> 

SW  S-wave  Velocity 
<km/sec> 

1 

2.10 

2.26 

2 

1.00 

3.50 

3 

1.00 

3.52 

4 

2.10 

3.57 

5 

2.00 

3.58 

6 

2.25 

3.61 

7 

2.25 

3.69 

8 

2.25 

3.70 

9 

2.00 

3.39 

10 

2.75 

3.51 

11 

2.25 

3.60 

12 

2.90 

3.79 

13 

1.40 

3.93 

14 

2.00 

3.94 

15 

2.50 

3.80 

16 

2.50 

3.54 

17 

2.80 

3.78 

18 

3.00 

4.20 

■-9 

3.00 

4.38 

20 

2.00 

4.18 

21 

2.00 

4.24 

12 

2.00 

4.22 

23 

2.00 

4.13 

24 

2.50 

4.02 

25 

2.50 

4.51 

26 

2.50 

4.46 

27 

2.50 

4.56 

28 

2.50 

4.51 

29 

0.00 

4.44 
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near  52  km  depth  which  produces  the  high  amplitude  phases 
between  22  to  25  seconds  in  the  synthetic  receiver  function, 
the  associated  errors  make  interpreting  the  abruptness  of 
this  boundary  difficult. 

GEOLOGIC  INTERPRETATION  OF  MODELING  RESULTS 

The  shear  wave  velocity  models  together  with  observed 
regional  compressional  velocities,  Rio  Grande  rift  geotherms 
(Decker  and  Smithson,  1975)  and  igneous  and  metamorphic  phase 
diagrams  (Figures  26  and  27)  comprised  the  basis  for  these 
lithological  determinations.  Absolute  lithology 
determination  is  supported  by  information  gathered  by  N.  I. 
Christensen  (1982)  on  characteristic  rock  velocities  and  by 
A.  W,  Hatheway  and  G.  A.  Kiersch  (1982)  on  engineering 
properties  of  rocks.  Upper  mantle  structural  interpretations 
are  made  here  in  light  of  the  observed  volcanism  in  the 
central  rift  region,  together  with  various  references  on 
experimentally  derived  upper  mantle  models  using  xenolith 
studies . 

UPPER  CRUST  i_<  15  KM  DEPTH)  ; 

Beneath  the  near  surface  sedimentary  layer,  the  upper 
crustal  basement  beneath  ANMO  is  known  to  be  comprised  of 
Precam, brian  granitic  rocks  (Vp  =  5.9  -  6.1  km/sec).  Shear 
velocities  for  all  vhree  back  azimuths  range  from  3.4  km/sec 
near  the  surface  to  3.5  -  3.7  km/sec  near  15  km  depth 
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(Poisson  ratios  ranging  between  .21  to  .24).  The  gradational 
character  of  the  shear  velocities  may  indicate  a  transition 
from  more  felsic  to  intermediate  composition  granitic  roc)cs. 

MIDDLE  CRUST  (15  -  25  ^  DEPTH) : 

Throughout  this  zone,  shear  velocities  in  the  models  are 
indicative  of  intermediate  composition  igneous  roc)cs  (e.g., 
granodiorite ) .  To  the  southeast,  velocities  grade  smoothly 
throughout  this  zone.  Higher  heat  flow  along  the  central 
rift  along  with  5  to  6  )cbars  pressure  at  15  to  20  km  depth 
would  tend  to  produce  a  slight,  reduction  in  the  shear 
velocity  of  approximately  2.7  gm/cc  density  granodiorite, 
which  we  observe  in  both  the  southwest  and  northwest  back 
azimuthal  models.  The  aforementioned  velocity  reversal  in 
the  southwest  back  azimuth  near  the  top  of  this  zone  could 
also  be  reasonably  explained  by  local  thermal  effects. 

Figure  28  shows  the  outline  of  the  Socorro  magma  body  as 
mapped  by  Rhinehart  et  al .  (1979).  The  propagation  path  for 
energy  arriving  from  the  Tonga  teleseisms  intersects  this 
horizon  (15  -  20  km  depth)  within  nearly  20  km  distance  from 
the  northernmost  extent  of  this  body.  However,  the  standard 
errors  in  velocity  mask  the  degree  to  which  these  thermal 
effects  may  attenuate  shear  velocities  in  this  vicinity.  The 
sout.hwest  model  demonstrates  a  relatively  steeper  velocity 
gradient  into  lower  crustal  depths  as  compared  to  the 
northwest  and  southeast  models. 
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LOWER  CRUST  AND  CRUST-MANTLE  TRANSITION  ( 25  -  40  ^  DEPTH ) ; 

In  each  model ,  an  increase  in  shear  velocity  is 
noted  within  the  upper  portion  of  this  zone.  This  may 
reflect  a  change  from  intermediate  igneous  lithologies  (e.g., 
granodiorite )  to  intermediate-to-maf ic  granulite  lithologies. 
These  anhydrous  granulites  can  form  in  continental  rift  zones 
under  the  influence  of  high  mantle-derived  heat  flow  (Weber, 
1986).  Figure  27  shows  that  for  the  Rio  Grande  rift  thermal 
and  pressure  regimes,  granulites  may  form  near  20  km  depth. 
According  to  Frost  and  Frost  (1987),  in  extensional  regimes, 
upwelling  basaltic  magmas  at  the  base  of  the  crust  introduce 
additional  heat  as  well  as  high  volumes  of  C02  which  act  to 
dehydrate  neighboring  rocks,  thereby  driving  granulite 
metamorphism  of  crustal  materials.  Along  the  rift  axis, 
shear  velocities  derived  here  of  approximately  3.75-3.80 
km/sec  and  compressional  velocities  near  6.5  km/sec  (Sinno  et 
al.,  1986;  Olsen,  1979)  are  consistent  with  2. 8-2. 9  gm/cc 
density  mafic  granulites.  In  addition,  seismic  reflection 
imaging  of  the  lower  crust  (e.g..  Brown  et  al.,  1979,  1980) 
has  revealed  layered,  discontinuous  reflectors  within  this 
zone  which  could  indicate  the  expected  layering  of  these 
granulites  with  mafic  rocks  and  preexisting  lower  velocity 
lithologies.  However,  the  layer  thickness  parameter izations 
of  2-3  km  used  in  this  present  study  fail  to  reveal 
lam.inations  in  the  1-D  models  locally  within  these  zones. 
Based  on  these  shear  velocities,  these  complex  zones  are 


Figure  26.  Southern  Rio  Grande  rift  geotherms :  qs  =  surface 
flux  assumed  in  temperature  calculations.  Wet 
solidus  curves:  G  =  granodior ite ,  P  =  peridotite, 
AB  =  alkaline  basalt,  OT  =  olivine  tholeiite. 

Dry  solidus  curves:  OQ  =  olivine  tholeiite  to 
quartz  eclogite  at  high  pressures,  OE  =  olivine 
tholeiite  to  olivine  eclogite  at  high  pressures, 
and  DP  ==  peridotite  (from  Decker  and  Smithson, 
1975  )  . 
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Figure  27.  Metamorphic  phase  diagram  with  Rio  Grande  rift 
geotherms  (from  Decker  and  Smithson,  1975) 
demonstrating  potential  for  lower  crustal 
granulite  facies  (modified  from  Ernst,  1976). 
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Proximity  of  (a)  propagation  path  of  seismic  energy 
arriving  from  Tonga  events  ( SW  BAZ)  to  (b)  mapped  extent 
of  Socorro  magma  body  (modified  from  Rhinehart  et  al . , 
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likely  interinediate-to-maf ic  granulites,  together  with 
underplated  mafics  and  other  precursor  crustal  rocks.  The 
refraction  work  of  Olsen,  et  al  (1979)  and  Sinno  et  al . 

(1986)  found  the  crust  mantle  boundary  between  31  and  34  km 
depth  along  the  southern  to  central  rift  axis.  This  zone 
corresponded  to  Pn  derived  velocities  between  7.6  -7.94 
km/sec,  with  7.77  km/sec  taken  as  the  average  Moho  v'elocity. 
However,  this  same  zone  corresponds  to  a  slight  attenuation 
in  the  shear  velocities.  Assuming  that  depth  related  errors 
are  minimal  between  these  analyses,  the  high  Poisson  ratios 
(near  .33)  evidenced  by  shear  wave  attenuation  within  this 
zone  could  be  attributed  to  the  presence  of  basaltic  material 
at  the  base  of  the  crust  from  partial  melting  within  the 
upper  mantle.  Plio-Pleistocene  volcanism  within  the  central 
rift  area  in  the  form  of  olivine  tholeiitic  basalts  indicate 
that  these  originated  from  partial  melting  within  the  upper 
mantle  (Baldridge,  1979).  This  mechanism  could  emplace  melt 
rich  in  silica,  olivine,  and  pyroxene  at  the  base  of  the 
crust,  which  may  produce  the  observed  shear  velocity  anomaly. 
However .  it  may  be  more  reasonable  to  assume  that  the  depth 
1.0  the  upper  mantle  in  the  vicinity  of  the  refraction 
profiles  would  differ  from  the  depths  to  the  upper  mantle 
noted  for  uhe  1-D  azimut.hal  models.  Then  this  zone  between 
30  and  34  km  in  the  models  with  attenuated  velocities  near 
3.6  to  3.7  km 'sec  would  not  correspond  l.o  the  refraction  Moho 
(l.e.,  would  lie  abov'e  it),  and  could  be  explained  as  an 


interval  characterized  by  partial  melting  of  the  lower  crust 
beneath  the  granulite  suite,  possibly  overlying  upwelled 
basalts  released  from  partial  melting  near  50  km  depth  in  the 
upper  mantle. 

Between  35  and  40  km  depth,  the  shear  velocities 
increase  rapidly  into  typical  upper  mantle  values  (4.25-4.35 
km/sec).  This  likely  reflects  a  transition  from  the 
overlying  lower  crustal  partial  melt  zone  to  the  high  silica- 
content  (near  50%)  olivine  and  pyroxene  upper  mantle  partial 
melt  (which  should  correspond  to  the  refraction  Moho,  with  Pn 
velocity  =  7.7  km, /sec )  to  spinel  peridotite  near  40  km  depth. 

UPPER  MANTLE  (BELOW  36-40  KM  DEPTH) ; 

The  velocities  within  the  upper  mantle  observed  in  the 
models  range  between  4.25  to  4.45  km/sec.  Given  the  higher 
rift  geotherm,  relative  to  adjacent  provinces,  these  are 
appropriate  velocities  for  upper  mantle  between  about  40  and 
60  km  depth,  where  a  peridotite  composition  with  some  partial 
melt  present  should  prevail.  Poisson  ratios  based  on  these 
determination  and  compress: onal  velocities  calculated  oy 
Sinno  (1984)  have  an  average  range  between  .28  and  ,31.  The 
pressure  gradient  within  this  zone  (approximately  .27 
kbar/km)  may  alone  account  for  the  slight  positive  shear 
velocity  gradient.  Although  the  velocit.y  step  near  52  km 
depth  obser\ned  in  the  southwest  model  does  produce  a 
satisfactory  waveform,  m,atch  for  phases  between  22  and  25 
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seconds  in  the  southwest  receiver  function,  one  cannot  be 
confident  that  these  are  not  intralayer  muluiple  phases  from 
the  crust-mantle  transition  which  cannot  be  generated  using 
the  fast  raytracing  method  of  Langston.  If  such  an  anomaly 
does  exist,  however,  partial  melt  between  45  and  50  km  could 
explain  this  shear  wave  attenuation.  This  depth  could  also 
mark  the  zone  where  replacement  of  spinel  by  garnet  in  upper 
mantle  peridot Ite  occurs,  although  it  is  j ^certain  if  the 
associated  density  change  would  be  sufficient  to  produce  this 
velocity  anomaly.  Another  possible  expxanation  could  be  that 
this  zone  marks  the  lithospher.  -asthenosphere  transition 
zone . 

EPT-ALQ  SURFACE  WAVE  DISPERSION  STUDY 

Table  6  lists  the  events  used  as  sources  in 
determination  of  int3rstation  dispersion.  These  earthquakes 
were  located  in  the  eastern  Pacific  between  1980  and  1982. 

The  same  EPT  and  ALQ  analog  recorded  long  period  vertical 
component  s'^ismograras  from  the  1986  study  were  obtained  and 
peaks,  troughs,  and  inflection  points  of  the  surface  wave 
trains  were  redigitized.  The  distance  from  the  event  to  the 
station,  relative  start  time  (0  seconds)  and  end  time  for 
digiti  zauion ,  and  time  in  seconds  from  the  e',en_  time  to  the 
start  of  the  record  were  recorded  for  velocity  determ.ination . 
The  raw  digitized  seismograms  were  interpolated  to  1 
samiole/sec  with  a  piecewise,  cubic,  polynomial  fit  (Wiggins, 
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TABLE  6:  Source  Events  and  EPT-ALQ  Station  Data 


Date 

Origin  Time 

Lat .  Lon . 

<degrees> 

Azimuth 
to  EPT 
<degrees> 

Dist . 
to  EPT 
<km> 

09/29/80 

2221:45.3 

18.33  N 

106.37  W 

-0.5 

1489.1 

05/05/82 

1200:41.7 

22.09  N 

107.24  W 

3.7 

1075.3 

12/20/81 

0907:32.3 

4.40  S 

105.30  W 

-1.8 

4001.0 

Station 

Lat. 

Lon. 

Azimuth 
to  EPT 
<degrees> 

Dist . 
to  EPT 
<km> 

EPT 

31.77  N 

106.51 

w  - 

- 

— 

ALQ 

34.94  N 

106.46 

W  1.0 

351.6 
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1976)  prior  to  filtering. 

Figures  29  through  34  demonstrate  filtering  of  the  EPT 
and  ALQ  seismograms.  In  each,  part  (a)  is  the  windowed 
seismogram.  A  quasi-pink  filter  (QPF)  was  applied  to  the 
Fourier  transform  of  each  signal  in  order  to  boost  long 
period  gain  lost  due  to  WWSSN  instrument  response  for  low 
frequencies  (Dean  and  Keller,  1990).  Part  (b)  shows  the 
multiple  filter  technique  (MFT)  derived  group  dispersion 
(triangles),  and  phase  match  iterated  (PMI)  phase  (solid)  and 
group  (dashed)  curves.  For  each  station  analysis,  the  time 
variable  filtered  (TVF)  seismograms  appear  in  part  (c). 

The  windowed  seismogram  in  Figures  29  (from  EPT)  and  30 
(from  ALQ)  are  for  an  event  off  the  coast  of  Jalisco,  Mexico 
on  September  29,  1980.  Both  signals  are  fairly  clean,  with 
the  exception  of  the  long  period  noise  component  with 
apparent  velocity  between  2.1  and  2.4  km/sec  in  the  ALQ 
recording.  PMI  results  in  essentially  no  refinement  in  group 
velocity  (Figures  29(b)  and  30(b)).  Figures  29(c)  and  30(c) 
are  the  signals  after  time  variable  filtering.  Later 
arriving  interference  (presumably  multipath)  was  effectively 
removed  in  the  TVF  process. 

Figures  31(a)  and  32(a)  are,  respectively,  the  windowed 
EPT  and  ALQ  seismograms  for  the  May  5,  1982  event  off  the 
coast  of  central  Mexico.  A  higher  mode  early  arrival 
(apparent  velocity  between  3.2  and  3.7  km/sec)  is  noted  in 
both  signals.  PMI  applications  for  the  EPT  signal  using  an 
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experimental  range  of  integration  constants  between  40  to  60 
seconds  resulted  failed  to  produce  any  refinement  to  MFT 
selected  group  velocities  between  40  and  50  second  periods 
(Figure  31(b)),  where  clearly  a  discrepancy  in  PMI  determined 
phase  and  group  velocity  exists.  In  both  TVF  seismograms, 
higher  mode  and  later  arriving  interference  are  significantly 
reduced  (Figures  31(c)  and  32(c)). 

In  Figures  33  and  34,  filtering  results  for  the  EPT  and 
ALQ  signals  from  a  December  20,  1989  event  in  the  northern 
Easter  Island  cordillera  event  are  shown.  PMI  results  in  a 
slight  refinement  in  MFT  selected  group  velocities.  Again, 
higher  mode  propagation  and  later  arriving  interference  has 
been  removed  in  the  TVF  process. 

Filtering  of  the  interstation  Green's  functions  for  the 
three  events  is  shown  in  Figures  35,  36,  and  37.  In  each, 
part  (a)  is  the  windowed  cross  correlation  of  the  EPT  and  ALQ 
filtered  signals.  Interstation  dispersion  is  shown  in 
part(b).  The  extreme  deviation  in  the  PMI  group  velocity 
from  the  MFT  group  velocity  noted  in  the  EPT-ALQ  dispersion 
for  the  Easter  Island  event  (Figure  37(b))  deserves  comment. 
The  MFT  selected  velocities  in  this  period  range  were  edited 
out  because  they  looked  incorrect.  As  a  result,  the  PMI 
application  yielded  this  refinement  to  the  group  velocity 
curve,  which  is  close  to  the  original  MFT  selection. 

Figures  38,  39,  and  40  illustrate  the  interstation 
dispersion  determined  from  the  filtered  Green's  functions. 
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MFT  group  velocities  are  shown  in  part  (a),  with  the  fast 
(solid)  and  slow  (dashed)  group  velocity  curves.  PMI  for  the 
fast  and  slow  group  velocities  yields  fast  and  slow  phase 
velocities,  which  are  meaned  to  produce  the  phase  velocities 
(circles)  in  part  (b).  The  curves  in  part  (b)  are  the 
simultaneous  least-squares  fit  to  the  mean  phase  and  MFT 
group  data.  The  self-consistent  phase  and  group  velocities 
with  standard  errors  are  given  in  part  (c). 

Having  determined  the  EPT-ALQ  path  dispersion  for  each 
individual  event  (Figure  41(a)),  these  are  weighted  by  the 
inverse  square  of  their  standard  errors,  and  the  mean  group 
(triangles)  and  phase  (circles)  interstation  dispersion  is 
determined,  as  shown  in  Figure  41(b). 

In  Figure  42(b),  comparison  of  the  polynomial  fitted 
group  and  phase  dispersion  curves  from  Sinno  and  Keller's 
1986  study  (symbols)  with  the  simultaneous  fits  derived  here 
(curves)  reveals  close  agreement  in  phase  velocity  over  the 
10  to  65  second  period  range.  Group  velocity  deviations  are 
noted  at  the  extreme  periods  of  our  analysis.  The  techniques 
applied  here  have  resulted  in  a  smoother  estimal.ed  group 
dispersion  curve  over  the  period  range  of  analysis.  However, 
given  the  short  interstation  distance  (351  km)  between  EPT 
and  ALQ,  no  additional  dispersion  information  could  be  gained 
for  periods  outside  this  range. 
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Figures  29  -  34.  (following  pages'  Filtering  results  for 

individual  events;  (a)  Windowed,  resampled 
seismogram.  (b)  Dispersion  of  (a). 
Triangles  =  MFT  selected  group  velocities, 
dashed  line  =  PMI  group  velocities,  and 
solid  curve  =  PMI  phase  velocities,  (c)  TVF 
seismogram  based  on  PMI  group  velocity 
curve  in  (b) . 
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Figure  31. 
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EPT :  12/20/81  event 
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Figures  35  -  37.  (following  pages)  Filtering  results  for 

EPT-ALQ  interstation  Green's  functions. 

(a)  Windowed  cross  correlation  of  EPT  and 
ALQ  TVF  seismograms  for  single  events. 

(b)  Interstation  dispersion  of  seismogram 
in  (a).  Triangles  =  MFT  group  velocities, 
dashed  curves  =  PMI  group  velocities,  and 
solid  curves  =  PMI  phase  velocities. 

(c)  TVF  Green's  function  based  on  PMI  group 
velocity  curve  shown  in  (b) . 
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Figure  35. 
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EPT-ALQ:  5/5182 
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Figure  36. 
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Figures  38-40.  (following  pages)  EPT-ALQ  interstalion 

dispersion  determination  for  the  three 
events.  (a)  Triangles  =  MFT-selected  group 
velocities.  Solid  curves  =  fast  velocities, 
dashed  curves  =  slow  velocities,  which  are 
determined  by  standard  errors  of  the  least 
squares  smoothed  group  velocities  and  inputs 
to  PMI .  (b)  Triangles  =  MFT  group 

velocities,  circles  =  meaned  fast  and  slow 
PMI  output,  curves  =  simultaneous  fit 
dispersion.  (c)  Final  phase  and  group 
dispersion  with  standard  errors. 
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EPT-ALQ:  (5/^162)  dispersion 
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EP  T-ALQ:  (12/20  /8/)  dispersion 


Figure  40. 
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EPT-ALQ:  dispersion 
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igure  41.  Simultaneous  fit  of  EPT-ALQ  interstation 

dispersion  for  the  three  events.  (a)  Individual 
event  dispersion.  (b)  Weighted  means  (triangles 
group  velocity,  circles  =  phase  velocity)  for  the 
three  events  with  error  bars.  Curves  are  the 
simultaneous  fit  of  the  m»-an  data. 
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EPT'ALQ:  final  dispersion 
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Figure  42.  Comparison  of  final  self -consistent,  interstation 
dispersion  with  results  from  the  1986  study. 

(a)  Final  interstation  dispersion,  this  study. 
Triangles  =  group  velocities,  circles  =  phase 
velocities,  with  error  bars.  (b)  Symbols  =  1986 
study  interstation  dispersion  results,  curves  = 
final  interstation  dispersion,  this  study. 
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MODELING  OF  GROUP  AND  PHASE  DISPERSION 

Implementing  the  first  order  difference  inversion,  two 
starting  models  are  used  in  the  determination  of  interstation 
shear  velocity  structure.  First  of  all,  a  constrained 
starting  model  based  on  the  final  model  from  the  southwest 
back  azimuth  receiver  function  analysis  is  constructed.  This 
biased  model  is  a  smoothed  version  of  the  southwest  model, 
confined  within  the  standard  error  bounds  in  the  derived 
velocity  model  (Figure  43).  Poisson  ratios  of  .24  for  the 
upper  crust,  .27  for  the  lower  crust,  and  .29  for  the  upper 
mantle  are  assigned.  Initial  P-wave  velocities  are  then 
determined  based  on  the  relationship  of  S-wave  velocity  and 
Poisson  ratio. 

The  derived  S-wave  velocity  model  (Figure  44,  and  Table 
7)  generates  the  dispersion  (solid  symbols)  shown  in  Figure 
45  and  listed  in  Table  8.  This  demonstrates  a  significant 
match  to  the  group  and  phase  dispersion  (open  symbols) 
determined  from  the  filtered  interstation  Green's  functions. 
The  model  of  shear  structure  derived  from  surface  wave 
dispersion  should  be  expected  to  deviate  slightly  from  the 
local  models  for  structure  beneath  ANMO  in  that  (a)  it 
represents  velocity  structure  averaged  over  the  interstation 
distance  through  the  rift  province,  and  (b)  dispersion 
results  are  most  reliable  for  the  mid-periods;  hence,  uppe  ■ 
and  lower  portions  of  the  derived  models  are  less  reliable. 
Average  shear  velocities  near  3.4  km/sec  for  the  upper  crust. 


ID  b:  R  D" 
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Figure  43.  EPT-ALQ  S-wave  velocity  starting  model  for 

surface  wave  inversion,  constrained  with  derived 
model  from  southwest  back  azimuth  receiver 
function  modeling.  Solid  line  =  starting  model, 
dashed  lines  =  +/-  2  standard  errors  in  southwest 
model . 


H  j_  cd  Id  a 


Inversion  results  from  starting  velocity  model 
outlined  in  Figure  43.  Solid  line  =  derived 
model,  dashed  lines  =  +/-  2  standard  errors  in 
soutnwest  model . 


Figure  44 . 
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TABLE  7:  Final  S-wave  Velocity  Model —  EPT-ALQ  Inter station 
Path  (Inversion  Constrained  with  SW  Back  Azimuth 
Receiver  Function  Results ) 


Layer  # 

Thickness 

<km> 

EPT-ALQ  S-wave 
<km/sec> 

1 

2.10 

2.09 

2 

1.00 

3.43 

3 

1.00 

3 .45 

4 

2.10 

3.32 

5 

2.00 

3.24 

6 

2.25 

3.23 

7 

2.25 

3.36 

8 

2.25 

3.49 

9 

2.00 

3.58 

10 

2.75 

3.75 

11 

2.25 

3.76 

12 

2.90 

3.95 

13 

1.40 

3.86 

14 

2.00 

3.85 

15 

2.50 

3.81 

16 

2.50 

3.86 

17 

2.80 

4.00 

18 

3.00 

4.10 

19 

3.00 

4.06 

20 

2.00 

4.11 

21 

2.00 

4.10 

22 

2.00 

4 . 10 

23 

2.00 

4 . 06 

24 

2.50 

3.94 

25 

2.50 

3.82 

26 

2.50 

4 .25 

27 

2.50 

4.27 

28 

0.00 

4  .  14 

100 


TABLE  8:  Final  Observed-Synthetic  Dispersion —  EPT-ALQ 

Interstation  Path  (Synthetic  Dispersion  for  Derived 
Model  from  Inversion  Constrained  with  SW  Back 
Azimuth  Receiver  Function  Results) 


Period 

<sec> 

Simultaneous  Fit 
<km/sec> 

Synthetic 

<km/sec> 

Phase 

Phase 

10.56 

3.060  +/- 

.024 

3.090  ■ 

12.13 

3.160  +/- 

.016 

3.166 

13.92 

3.250  +/- 

.012 

3.245 

16.00 

3.331  +/- 

.011 

3.324 

18.38 

3.402  +/- 

.012 

3.397 

21.11 

3.464  +/- 

.013 

3.464 

24.25 

3.519  +/- 

.015 

3.522 

27.89 

3.566  +/- 

.018 

3.573 

32.00 

3.608  +/- 

.021 

3.615 

36.76 

3.647  +/- 

.025 

3.651 

42.22 

3.684  +/- 

.030 

3.680 

48 . 50 

3.721  +/- 

.037 

3.704 

55.72 

3,762  +/- 

.049 

3.724 

64.00 

3.808  +/- 

.069 

3.741 

Group 

Group 

10.56 

2.458  +/- 

.072 

2.642 

12.13 

2.597  +/- 

.054 

2.685 

13.93 

2.731  +/- 

.038 

2.756 

16 . 00 

2.860  +/- 

.028 

2.849 

18.38 

2.980  +/- 

.024 

2.955 

21.11 

3.090  +/- 

.025 

3.064 

24.25 

3.188  +/- 

.026 

3.170 

27 . 86 

3.271  +/- 

.025 

3.267 

32.00 

3.341  +/- 

.027 

3.355 

36.76 

3.395  +/- 

.035 

3.432 

42.22 

3.436  +/- 

.051 

3.498 

48 . 50 

3.462  +/- 

.072 

3.552 

55.72 

3.477  +/- 

.097 

3.597 

64.00 

3.481  +/- 

.123 

3.634 

101 


3.6  km/sec  for  the  mid-crust,  and  3. 7-3. 8  for  the  lower  crust 
are  seen  in  this  model.  In  light  of  the  above  considerations, 
these  results  agree  satisfactorily  with  the  southwest  back 
azimuth  modeling  results. 

A  second  approach  at  resolving  interstation  shear 
structure  involved  beginning  with  a  three  layer  model  with 
2.25  km/sec  for  the  near  surface,  3.4  km/sec  for  -the  entire 
crust,  and  4.0  km/sec  for  the  upper  mantle  (Figure  46). 
Inversion  of  this  unbiased  starting  model  resulted  in  a 
convergence  toward  the  previously  derived  models  (Figure  47). 
Although  this  derived  model  contains  unrealistic  components 
(for  example,  the  velocity  reversal  at  4  to  10  km  depth,  and 
the  strong  first  order  crust-mantle  boundary  imposed  in  the 
starting  model),  one  can  observe  that  the  mid-to-lower  crust 
portion  of  this  model  moves  within  the  bounds  of  the  model 
determined  from  receiver  function  analysis.  The 
corresponding  dispersion  shown  in  Figure  48  demonstrates  the 
significant  fit  to  observed  dispersion.  This  further 
supports  the  previously  presented  1-D  determinations  of  shear 
wave  velocity  ,  structure  beneath  ANMO  in  the  Rio  Grande  rift. 


CONCLUSION 

The  ANMO  receiver  function  study  results  are  significant 
in  at  least  two  ways.  Firstly,  it  has  been  shown  here  that 
conventional  methods  of  modeling  detailed  receiver  structure 
can  be  applied  for  three  component  digital  seismic  stations 


X 


Figure  46.  3-layer  EPT-ALQ  S-wave  velocity  starting  model 
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S  —  VEILOCTTY  <k.m^sec::> 


2.0  3.0  ‘4.0  5.0 


Figure  47.  Inversion  results  from  starting  model  outlined  in 
Figure  46,  demonstrating  convergence  toward 
southwest  back  azimuth  model  in  the  mid-to-low^r 
crustal  range. 
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lacking  broadband  instrumentation.  Merging  the  long  and 
short  period  digital  records  provides  the  appropriate 
bandwidth  for  such  high  resolution  analysis.  Secondly,  more 
information  on  deep  shear  structure  within  the  central  Rio 
Grande  rift  has  been  provided.  Integrating  these  local 
findings  with  results  from  future  such  studies  within  the  Rio 
Grande  rift  (e.g.,  EPT  receiver  structure  in  the  southern 
rift,  and  local  studies  using  portable  broadband  recording 
instruments)  should  aid  in  detailed  regional  mapping  of  the 
extremely  heterogeneous  lower  crustal  and  upper  mantle 
structure  within  this  province. 

The  EPT-ALQ  dispersion  analysis  results  are  quite 
consistent  with  those  of  the  original  study.  This  is 
understandable,  given  that  information  outside  the  original 
period  range  analyzed  could  not  be  extracted  due  to  the  band- 
limitedness  of  the  dispersion  data  for  such  a  short 
interstation  path.  However,  this  study  provided  another 
estimate  of  average  shear  velocity  structure  along  the  rift 
axis  using  the  same  data,  but  processed  with  different 
filtering  techniques.  Also,  these  S-wave  velocity  estimates 
provide  a  satisfactory  tie  with  the  mid-to-lower  crustal 
portion  of  the  velocity  model  derived  in  the  southwest  back 
azimuth  receiver  function  analysis.  These  results,  together 
wiLh  those  from  previous  Rio  Grande  rift  geophysical  studies, 
will  hopefully  provide  both  reliable  control  for  future 
seismic  modeling,  and  a  better  understanding  of  the 
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petrologic  processes  taking  place  within  the  rift. 


APPENDIX 


LIST  OF  COMPUTER  PROGRAMS 
USED  IN  THESIS. 
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PROGRAM:  RAY3D 

AUTHOR:  Dr.  T.J.  Owens 

MODIFICATIONS  BY  B.P.  MURPHY:  Replaced  S.A.C.  (Seismic 

Analysis  Code)  input/output 
subroutines  with  HP-9000 
compatible  routines,  which 
manipulate  S.A.C.  format  ASCII 
files . 

PROGRAM  PURPOSE:  Implementation  of  Langston's  3-D  raytracing 

method,  used  in  modeling  receiver 
function  waveforms. 


PROGRAM:  PWAVEQN 

AUTHOR:  Dr.  T.J.  Owens 

MODIFICATIONS  BY  B.P.  MURPHY:  Replaced  S.A.C.  input/output 

subroutine  calls. 

PROGRAM  PURPOSE:  Implementation  of  Langston's  source 

equalization  deconvolution  method. 
Deconvolves  vertical  component  from 
horizontal  component  (radial  and 
tangential)  seismograms  to  produce  receiv’^er 
functions . 


PROGRAM:  STACK 

AUTHOR:  Dr.  T.J.  Owens 

MODIFICATIONS  BY  B.P.  MURPHY:  Added  code  for  determining 

standard  deviations  (+/-  1 
S.D.  about  the  mean  stac}ced 
trace)  based  on  data  variance. 
Upper  and  lower  bound  output 
traces  are  produced. 

PROGRAM  PURPOSE:  Determines  a  mean  receiver  function  trace 

for  a  particular  bac)c  azimuth  by  stac)cing  a 
variable  number  of  receiver  functions  with 
approximately  equal  source  distances  and 
directions . 
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PROGRAM;  INVINI 

AUTHOR:  Dr.  T.J.  Owens 

MODIFICATIONS  BY  B.P.  MURPHY:  Replaced  S.A.C.  input/output 

subroutine  calls. 

PROGRAM  PURPOSE:  Initializes  a  binary  output  file,  which  is 

input  file  for  the  inversion  program 
(TIMINV).  This  file  contains  all  user 
selected  inversion  options,  plus  user 
specified  starting  model  and  receiver 
function  trace. 


PROGRAM:  TIMINV 

AUTHOR:  Dr.  T.J.  Owens 

MODIFICATIONS  BY  B.P.  MURPHY:  Replaced  S.A.C.  input/output 

subroutine  calls. 

PROGRAM  PURPOSE:  Given  the  binary  parameter  output  file  from 

INVINI,  TIMINV  inverts  a  receiver  function 
trace  for  S-wave  velocity.  RAY3D  is  called 
for  forward  calculations.  Updated  S-wave 
velocity  model  with  upper  and  lower  bound 
velocity  models  and  the  corresponding 
synthetic  radial  component  receiver 
function  are  the  output  produced. 


PROGRAM;  HPMERGE 

AUTHOR;  B.  Belmont 

MODIFICATIONS  BY  B.  P.  MURPHY:  None 

PROGRAM  PURPOSE;  Implemetation  of  Harvey  and  Choy ' s  data 

merging  technique.  Merges  WWSSN  digital 
three  component  LP  (long  period)  and  SP 
(short  period)  seismograms  to  produce  BB 
(broadband)  seismograms  for  receiver 
function  waveform  modeling. 


PROGRAM:  ROTATE  (subroutine) 

AUTHOR:  Dr.  D.  I.  Doser 

MODIFICATIONS  BY  B.  P.  MURPHY:  Coded  into  main  calling 

routine,  ROTE. 

PROGRAM  PURPOSE;  Rotates  N-S  and  E-W  seismograms  in  source 

bac)^  azimuth  to  produce  radial  and 
tangential  component  seismograms. 
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PROGRAM:  SEISMO 

AUTHOR:  Dr.  E.  A.  Dean 

MODIFICATIONS  BY  B.  P.  MURPHY:  NONE 

PROGRAM  PURPOSE:  Interactive  processina  routine  for 

determining  group  and  phase  velocities  and 
their  standard  errors  for  multi-event 
surface  wave  data. 


PROGRAM:  VS 3 

AUTHOR:  B.  P.  MURPHY 

PROGRAM  PURPOSE:  Plots  S-wave  velocity  vs.  Depth  models 

for  receiver  function  analysis.  Standard 
error  bars  may  be  included  in  the  plot. 
VELPARAMS  is  the  required  control  file. 


PROGRAM:  S_GRAPH 

AUTHOR:  B.  P,  MURPHY 

PROGRAM  PURPOSE:  Plots  a  user  specified  number  of 

seismograms  stored  in  S.A.C.  ASCII  format. 
S_CONFIG  is  the  required  control  file. 


PROGRAM:  FIGURE 

AUTHORS:  D.  G.  Roberts  and  Dr.  E.  A.  Dean 

MODIFICATIONS  BY  B.  P.  MURPHY:  None 

PROGRAM  PURPOSE:  Plots  SEISMO  output  (dispersion  curves  and 

filtered  seismograms). 


PROGRAM:  SURF 

AUTHOR:  Dr.  D.  R.  Russell 

MODIFICATIONS  BY  B.  P.  MURPHY:  None 

PROGRAM  PURPOSE:  Surface  wave  inversion  program.  Output  is 

updated  S-wave  velocity  model,  and 
synthetic  group  and  phase  dispersion  data. 
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